/*-*-mode: 

!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/* fun-pw91x.c:

   NOTE: This is the equivalent of the pwggaIIx functional.

   Automatically generated code implementing PW91X functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional PW91xFunctional;" to 'functionals.h'
    2. add "&PW91xFunctional," to 'functionals.c'
    3. add "fun-PW91x.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
rho:  rhoa + rhob;
grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
zeta: (rhoa-rhob)/(rhoa+rhob);

kfa: (6*%PI*%PI*rhoa)^(1/3);
kfb: (6*%PI*%PI*rhob)^(1/3);
exa: -3*kfa/(4*%PI);
exb: -3*kfb/(4*%PI);

sa:  grada/(2*kfa*rhoa);
sb:  gradb/(2*kfb*rhob);
F0a: 1+0.19645*sa*asinh(7.7956*sa);
F0b: 1+0.19645*sb*asinh(7.7956*sb);
F1a: F0a + (0.2743-0.1508*exp(-100*sa^2))*sa^2;
F1b: F0b + (0.2743-0.1508*exp(-100*sb^2))*sb^2;
F2a: F0a + 0.004*sa^4;
F2b: F0b + 0.004*sb^4;
Fa: F1a/F2a;
Fb: F1b/F2b;

Exa: 2*rhoa*exa*Fa;
Exb: 2*rhob*exb*Fb;

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb)/2;


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"

/* INTERFACE PART */
static integer pw91x_isgga(void) { return 1; }
static integer pw91x_read(const char *conf_line);
static real pw91x_energy(const FunDensProp* dp);
static void pw91x_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void pw91x_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void pw91x_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void pw91x_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional PW91xFunctional = {
  "PW91x",       /* name */
  pw91x_isgga,   /* gga-corrected */
   1,
  pw91x_read,
  NULL,
  pw91x_energy,
  pw91x_first,
  pw91x_second,
  pw91x_third,
  pw91x_fourth
};
 
/* IMPLEMENTATION PART */
static integer
pw91x_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
pw91x_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/t1;
    t4 = 1/pow(3.141592653589793,2.666666666666667);
    t5 = 1/pow(3.141592653589793,0.666666666666667);
    t6 = pow(rhoa,1.333333333333333);
    t7 = 1/t6;
    t8 = 0.098225*t3*t5*grada*asinh(3.8978*t3*t5*grada*t7)*
        t7;
    t9 = 1/pow(6.0,0.666666666666667);
    t10 = 1/pow(3.141592653589793,1.333333333333333);
    t11 = pow(grada,2.0);
    t12 = 1/pow(rhoa,2.666666666666667);
    t13 = pow(rhob,1.333333333333333);
    t14 = 1/t13;
    t15 = 0.098225*t3*t5*gradb*asinh(3.8978*t3*t5*gradb*t14)*
        t14;
    t16 = pow(gradb,2.0);
    t17 = 1/pow(rhob,2.666666666666667);

   /* code */
    res = 0.5*(-1.5*t1*t13*t2*(0.25*(0.2743-0.1508/pow(2.718281828459045,
        25.0*t10*t16*t17*t9))*t10*t16*t17*t9+t15+1.0)/(4.166666666666667E-5*
        t3*t4*pow(gradb,4.0)/pow(rhob,5.333333333333333)+t15+1.0)-
        1.5*t1*t2*t6*(0.25*(0.2743-0.1508/pow(2.718281828459045,25.0*
        t10*t11*t12*t9))*t10*t11*t12*t9+t8+1.0)/(4.166666666666667E-5*
        t3*t4*pow(grada,4.0)/pow(rhoa,5.333333333333333)+t8+1.0));

    return res;
}

static void
pw91x_first_helper(real rhoa, real grada, real *res)
{    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/t1;
    t4 = 1/pow(3.141592653589793,2.666666666666667);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/pow(3.141592653589793,0.666666666666667);
    t8 = pow(rhoa,1.333333333333333);
    t9 = 1/t8;
    t10 = asinh(3.8978*t3*t7*grada*t9);
    t11 = 0.098225*t3*t7*grada*t10*t9;
    t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
    t13 = 1/t12;
    t14 = 1/pow(6.0,0.666666666666667);
    t15 = 1/pow(3.141592653589793,1.333333333333333);
    t16 = pow(grada,2.0);
    t17 = 1/pow(rhoa,2.666666666666667);
    t18 = 1/pow(2.718281828459045,25.0*t14*t15*t16*t17);
    t19 = 0.2743-
        0.1508*t18;
    t20 = 0.25*t14*t15*t16*t17*t19+t11+1.0;
    t21 = 1/pow(rhoa,6.333333333333333);
    t22 = 1/sqrt(15.19284484*t14*t15*t16*t17+1.0);
    t23 = 1/pow(rhoa,3.666666666666667);
    t24 = -0.510481873333333*t14*t15*t16*t22*t23;
    t25 = -0.130966666666667*t10*t3*t7*grada/pow(rhoa,2.333333333333333);
    t26 = 1/
        pow(t12,2.0);
    t27 = pow(grada,3.0);
    t28 = 0.382861405*t14*t15*grada*t22*t17;
    t29 = 0.098225*t3*t7*t10*t9;

   /* code */
    res[0] = 0.5*(-2.0*t1*t13*t2*t20*pow(rhoa,0.333333333333333)+
        1.5*t1*t2*t20*(t25+t24-2.222222222222222E-4*t3*t4*t5*t21)*
        t26*t8-1.5*t1*t13*t2*(t25+t24-0.666666666666667*t14*t15*t16*
        t19*t23-0.418888888888889*t3*t4*t5*t21*t18)*t8);
    res[1] = 0.5*(1.5*t1*t2*t20*t26*(t29+t28+1.6666666666666666E-4*
        t3*t4*t27*t6)*t8-1.5*t1*t13*t2*t8*(0.5*t14*t15*t17*t19*grada+
        t29+t28+0.314166666666667*t3*t4*t27*t6*t18));
}

static void
pw91x_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    pw91x_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
pw91x_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/t1;
    t4 = 1/pow(3.141592653589793,2.666666666666667);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/pow(3.141592653589793,0.666666666666667);
    t8 = pow(rhoa,1.333333333333333);
    t9 = 1/t8;
    t10 = asinh(3.8978*t3*t7*grada*t9);
    t11 = 0.098225*t3*t7*grada*t10*t9;
    t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
    t13 = 1/t12;
    t14 = pow(rhoa,0.333333333333333);
    t15 = 1/pow(6.0,0.666666666666667);
    t16 = 1/pow(3.141592653589793,1.333333333333333);
    t17 = pow(grada,2.0);
    t18 = 1/pow(rhoa,2.666666666666667);
    t19 = 1/pow(2.718281828459045,25.0*t15*t16*t17*t18);
    t20 = 0.2743-
        0.1508*t19;
    t21 = 0.25*t15*t16*t17*t18*t20+t11+1.0;
    t22 = 1/pow(rhoa,6.333333333333333);
    t23 = sqrt(15.19284484*t15*t16*t17*t18+1.0);
    t24 = 1/t23;
    t25 = 1/pow(rhoa,3.666666666666667);
    t26 = -0.510481873333333*t15*t16*t17*t24*t25;
    t27 = 1/pow(rhoa,2.333333333333333);
    t28 = -0.130966666666667*t3*t7*grada*t10*t27;
    t29 = t28+t26-2.222222222222222E-4*t3*t4*t5*t22;
    t30 = 1/pow(t12,2.0);
    t31 = t28+t26-0.666666666666667*t15*t16*t17*t20*t25-0.418888888888889*
        t3*t4*t5*t22*t19;
    t32 = pow(grada,3.0);
    t33 = 0.382861405*t15*t16*grada*t24*t18;
    t34 = 0.098225*t3*t7*t10*t9;
    t35 = t34+t33+1.6666666666666666E-4*t3*t4*t32*t6;
    t36 = 0.5*t15*t16*t18*t20*grada+t34+t33+0.314166666666667*
        t3*t4*t32*t6*t19;
    t37 = 1/pow(t12,3.0);
    t38 = 1/pow(rhoa,7.333333333333333);
    t39 = 1/pow(t23,3.0);
    t40 = -1.723482643374637*t3*t4*t5*t39*t38;
    t41 = 1/pow(rhoa,4.666666666666667);
    t42 = 2.552409366666666*t15*t16*t17*t24*t41;
    t43 = 0.305588888888889*t10*t3*t7*grada/pow(rhoa,3.333333333333333);
    t44 = 1/
        pow(3.141592653589793,4.0);
    t45 = 1.292611982530978*t3*t4*t32*t39*t22;
    t46 = -1.53144562*t15*t16*grada*t24*t25;
    t47 = -0.130966666666667*t3*t7*t10*t27;
    t48 = -0.969458986898234*t3*t4*t17*t39*t6;
    t49 = 0.76572281*t15*t16*t24*t18;

   /* code */
    res[0] = 0.5*(-1.5*t1*t13*t2*t31*t8+1.5*t1*t2*t21*t29*
        t30*t8-2.0*t1*t13*t14*t2*t21);
    res[1] = 0.5*(1.5*t1*t2*t21*t30*t35*t8-1.5*t1*t13*t2*
        t36*t8);
    res[2] = 0.5*(-1.5*t1*t13*t2*t8*(-4.65432098765432*t19*
        t44*pow(grada,6.0)/pow(rhoa,10.0)+t43+t42+2.444444444444445*
        t15*t16*t17*t20*t41+t40+3.77*t3*t4*t5*t38*t19)-0.666666666666667*
        t1*t13*t2*t21/pow(rhoa,0.666666666666667)+1.5*t1*t2*t21*t30*
        (t43+t42+t40+0.001407407407407*t3*t4*t5*t38)*t8-3.0*t1*t2*
        t21*pow(t29,2.0)*t37*t8+3.0*t1*t2*t29*t30*t31*t8-4.0*t1*t13*
        t14*t2*t31+4.0*t1*t14*t2*t21*t29*t30);
    res[3] = 0.5*(-1.5*t1*t13*t2*t8*(3.490740740740741*t19*
        t44*pow(grada,5.0)/pow(rhoa,9.0)-1.333333333333333*t15*t16*
        t20*t25*grada+t47+t46+t45-2.513333333333333*t3*t4*t32*t22*
        t19)+1.5*t1*t2*t21*t30*(t47+t46+t45-8.888888888888888E-4*t3*
        t4*t32*t22)*t8-3.0*t1*t2*t21*t29*t35*t37*t8+1.5*t1*t2*t29*
        t30*t36*t8+1.5*t1*t2*t30*t31*t35*t8-2.0*t1*t13*t14*t2*t36+
        2.0*t1*t14*t2*t21*t30*t35);
    res[4] = 0.5*(-1.5*t1*t13*t2*t8*(-2.618055555555555*t19*
        t44*t5/pow(rhoa,8.0)+t49+t48+0.5*t15*t16*t18*t20+1.570833333333333*
        t3*t4*t17*t6*t19)+1.5*t1*t2*t21*t30*(t49+t48+5.0E-4*t3*t4*
        t17*t6)*t8-3.0*t1*t2*t21*pow(t35,2.0)*t37*t8+3.0*t1*t2*t30*
        t35*t36*t8);

}

static void
pw91x_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    pw91x_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
pw91x_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63, t64, t65, t66;
    real t67, t68, t69, t70, t71, t72, t73, t74;
    real t75, t76, t77, t78, t79, t80, t81, t82;
    real t83;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/t1;
    t4 = 1/pow(3.141592653589793,2.666666666666667);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/pow(3.141592653589793,0.666666666666667);
    t8 = pow(rhoa,1.333333333333333);
    t9 = 1/t8;
    t10 = asinh(3.8978*t3*t7*grada*t9);
    t11 = 0.098225*t3*t7*grada*t10*t9;
    t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
    t13 = 1/t12;
    t14 = pow(rhoa,0.333333333333333);
    t15 = 1/pow(6.0,0.666666666666667);
    t16 = 1/pow(3.141592653589793,1.333333333333333);
    t17 = pow(grada,2.0);
    t18 = 1/pow(rhoa,2.666666666666667);
    t19 = 1/pow(2.718281828459045,25.0*t15*t16*t17*t18);
    t20 = 0.2743-
        0.1508*t19;
    t21 = 0.25*t15*t16*t17*t18*t20+t11+1.0;
    t22 = 1/pow(rhoa,6.333333333333333);
    t23 = sqrt(15.19284484*t15*t16*t17*t18+1.0);
    t24 = 1/t23;
    t25 = 1/pow(rhoa,3.666666666666667);
    t26 = -0.510481873333333*t15*t16*t17*t24*t25;
    t27 = 1/pow(rhoa,2.333333333333333);
    t28 = -0.130966666666667*t3*t7*grada*t10*t27;
    t29 = t28+t26-2.222222222222222E-4*t3*t4*t5*t22;
    t30 = 1/pow(t12,2.0);
    t31 = t28+t26-0.666666666666667*t15*t16*t17*t20*t25-0.418888888888889*
        t3*t4*t5*t22*t19;
    t32 = pow(grada,3.0);
    t33 = 0.382861405*t15*t16*grada*t24*t18;
    t34 = 0.098225*t3*t7*t10*t9;
    t35 = t34+t33+1.6666666666666666E-4*t3*t4*t32*t6;
    t36 = 0.5*t15*t16*t18*t20*grada+t34+t33+0.314166666666667*
        t3*t4*t32*t6*t19;
    t37 = 1/pow(rhoa,0.666666666666667);
    t38 = pow(t29,2.0);
    t39 = 1/pow(t12,3.0);
    t40 = 1/pow(rhoa,7.333333333333333);
    t41 = 1/pow(t23,3.0);
    t42 = -1.723482643374637*t3*t4*t5*t41*t40;
    t43 = 1/pow(rhoa,4.666666666666667);
    t44 = 2.552409366666666*t15*t16*t17*t24*t43;
    t45 = 1/pow(rhoa,3.333333333333333);
    t46 = 0.305588888888889*t3*t7*grada*t10*t45;
    t47 = t46+t44+t42+0.001407407407407*t3*t4*t5*t40;
    t48 = 1/pow(3.141592653589793,4.0);
    t49 = pow(grada,6.0);
    t50 = 1/pow(rhoa,10.0);
    t51 = t46+t44+2.444444444444445*t15*t16*t17*t20*t43+t42-
        4.65432098765432*t48*t49*t50*t19+3.77*t3*t4*t5*t40*t19;
    t52 = 1.292611982530978*t3*t4*t32*t41*t22;
    t53 = -1.53144562*t15*t16*grada*t24*t25;
    t54 = -0.130966666666667*t3*t7*t10*t27;
    t55 = t54+t53+t52-8.888888888888888E-4*t3*t4*t32*t22;
    t56 = pow(grada,
        5.0);
    t57 = 1/pow(rhoa,9.0);
    t58 = -1.333333333333333*t15*t16*t20*t25*grada+t54+t53+
        t52+3.490740740740741*t48*t56*t57*t19-2.513333333333333*t3*
        t4*t32*t22*t19;
    t59 = pow(t35,2.0);
    t60 = -0.969458986898234*t3*t4*t17*t41*t6;
    t61 = 0.76572281*t15*t16*t24*t18;
    t62 = t61+t60+5.0E-4*t3*t4*t17*t6;
    t63 = 1/pow(rhoa,8.0);
    t64 = t61+t60+0.5*t15*t16*t18*t20-2.618055555555555*t48*
        t5*t63*t19+1.570833333333333*t3*t4*t17*t6*t19;
    t65 = 1/pow(t12,4.0);
    t66 = 1/pow(t23,5.0);
    t67 = 1/pow(rhoa,11.0);
    t68 = -17.45640292348261*t48*t49*t66*t67;
    t69 = 1/pow(rhoa,8.333333333333334);
    t70 = 21.25628593495385*t3*t4*t5*t41*t69;
    t71 = 1/pow(rhoa,5.666666666666667);
    t72 = -13.49940953925926*t15*t16*t17*t24*t71;
    t73 = -1.01862962962963*t10*t3*t7*grada/pow(rhoa,4.333333333333333);
    t74 = 1/
        pow(3.141592653589793,5.333333333333333);
    t75 = 13.09230219261196*t48*t56*t66*t50;
    t76 = -13.35699048615344*t3*t4*t32*t41*t40;
    t77 = 6.295943104444444*t15*t16*grada*t24*t43;
    t78 = 0.305588888888889*t3*t7*t10*t45;
    t79 = -9.81922664445897*t48*t5*t66*t57;
    t80 = 7.755671895185867*t3*t4*t17*t41*t22;
    t81 = -2.041927493333334*t15*t16*t24*t25;
    t82 = 7.364419983344228*t48*t32*t66*t63;
    t83 = -3.877835947592934*t3*t4*grada*t41*t6;

   /* code */
    res[0] = 0.5*(-1.5*t1*t13*t2*t31*t8+1.5*t1*t2*t21*t29*
        t30*t8-2.0*t1*t13*t14*t2*t21);
    res[1] = 0.5*(1.5*t1*t2*t21*t30*t35*t8-1.5*t1*t13*t2*
        t36*t8);
    res[2] = 0.5*(-1.5*t1*t13*t2*t51*t8+1.5*t1*t2*t21*t30*
        t47*t8-3.0*t1*t2*t21*t38*t39*t8+3.0*t1*t2*t29*t30*t31*t8-0.666666666666667*
        t1*t13*t2*t21*t37-4.0*t1*t13*t14*t2*t31+4.0*t1*t14*t2*t21*
        t29*t30);
    res[3] = 0.5*(-1.5*t1*t13*t2*t58*t8+1.5*t1*t2*t21*t30*
        t55*t8-3.0*t1*t2*t21*t29*t35*t39*t8+1.5*t1*t2*t29*t30*t36*
        t8+1.5*t1*t2*t30*t31*t35*t8-2.0*t1*t13*t14*t2*t36+2.0*t1*t14*
        t2*t21*t30*t35);
    res[4] = 0.5*(-1.5*t1*t13*t2*t64*t8+1.5*t1*t2*t21*t30*
        t62*t8-3.0*t1*t2*t21*t39*t59*t8+3.0*t1*t2*t30*t35*t36*t8);
    res[5] = 0.5*(-1.5*t1*t13*t2*t8*(-310.2880658436214*t15*
        t19*t74*pow(grada,8.0)/pow(rhoa,13.66666666666667)+t73+t72-
        11.40740740740741*t15*t16*t17*t20*t71+t70+t68-31.74246913580246*
        t3*t4*t5*t69*t19+88.43209876543209*t48*t49*t67*t19)+0.444444444444444*
        t1*t13*t2*t21/pow(rhoa,1.666666666666667)+1.5*t1*t2*t21*t30*
        (t73+t72+t70-0.010320987654321*t3*t4*t5*t69+t68)*t8+9.0*t1*
        t2*t21*pow(t29,3.0)*t65*t8+4.5*t1*t2*t29*t30*t51*t8-9.0*t1*
        t2*t21*t29*t39*t47*t8+4.5*t1*t2*t30*t31*t47*t8-9.0*t1*t2*t31*
        t38*t39*t8-6.0*t1*t13*t14*t2*t51+6.0*t1*t14*t2*t21*t30*t47-
        12.0*t1*t14*t2*t21*t38*t39-2.0*t1*t13*t2*t31*t37+2.0*t1*t2*
        t21*t29*t30*t37+12.0*t1*t14*t2*t29*t30*t31);
    res[6] = 0.5*(-1.5*t1*t13*t2*t8*(232.716049382716*t15*
        t19*t74*pow(grada,7.0)/pow(rhoa,12.66666666666667)+4.88888888888889*
        t15*t16*t20*t43*grada+t78+t77+t76+t75-59.34259259259259*t48*
        t56*t50*t19+18.15185185185185*t3*t4*t32*t40*t19)+1.5*t1*t2*
        t21*t30*(t78+t77+t76+0.00562962962963*t3*t4*t32*t40+t75)*t8+
        9.0*t1*t2*t21*t35*t38*t65*t8+3.0*t1*t2*t29*t30*t58*t8-6.0*
        t1*t2*t21*t29*t39*t55*t8+3.0*t1*t2*t30*t31*t55*t8+1.5*t1*t2*
        t30*t35*t51*t8-3.0*t1*t2*t21*t35*t39*t47*t8+1.5*t1*t2*t30*
        t36*t47*t8-3.0*t1*t2*t36*t38*t39*t8-6.0*t1*t2*t29*t31*t35*
        t39*t8-4.0*t1*t13*t14*t2*t58+4.0*t1*t14*t2*t21*t30*t55-8.0*
        t1*t14*t2*t21*t29*t35*t39-0.666666666666667*t1*t13*t2*t36*
        t37+0.666666666666667*t1*t2*t21*t30*t35*t37+4.0*t1*t14*t2*
        t29*t30*t36+4.0*t1*t14*t2*t30*t31*t35);
    res[7] = 0.5*(-1.5*t1*t13*t2*t8*(-174.537037037037*t15*
        t19*t49*t74/pow(rhoa,11.66666666666667)+t81+t80+t79-1.333333333333333*
        t15*t16*t20*t25+38.39814814814815*t48*t5*t57*t19-9.215555555555554*
        t3*t4*t17*t22*t19)+1.5*t1*t2*t21*t30*t8*(t81+t80-0.002666666666667*
        t3*t4*t17*t22+t79)+9.0*t1*t2*t21*t29*t59*t65*t8+1.5*t1*t2*
        t29*t30*t64*t8-3.0*t1*t2*t21*t29*t39*t62*t8+1.5*t1*t2*t30*
        t31*t62*t8-3.0*t1*t2*t31*t39*t59*t8+3.0*t1*t2*t30*t35*t58*
        t8-6.0*t1*t2*t21*t35*t39*t55*t8+3.0*t1*t2*t30*t36*t55*t8-6.0*
        t1*t2*t29*t35*t36*t39*t8-2.0*t1*t13*t14*t2*t64+2.0*t1*t14*
        t2*t21*t30*t62-4.0*t1*t14*t2*t21*t39*t59+4.0*t1*t14*t2*t30*
        t35*t36);
    res[8] = 0.5*(-1.5*t1*t13*t2*t8*(130.9027777777778*t15*
        t19*t56*t74/pow(rhoa,10.66666666666667)+t83+t82-23.5625*t48*
        t32*t63*t19+3.769999999999999*t3*t4*grada*t6*t19)+1.5*t1*t2*
        t21*t30*t8*(t83+0.001*t3*t4*grada*t6+t82)+9.0*t1*t2*t21*pow(t35,
        3.0)*t65*t8+4.5*t1*t2*t30*t35*t64*t8-9.0*t1*t2*t21*t35*t39*
        t62*t8+4.5*t1*t2*t30*t36*t62*t8-9.0*t1*t2*t36*t39*t59*t8);

}

static void
pw91x_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    pw91x_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
pw91x_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63, t64, t65, t66;
    real t67, t68, t69, t70, t71, t72, t73, t74;
    real t75, t76, t77, t78, t79, t80, t81, t82;
    real t83, t84, t85, t86, t87, t88, t89, t90;
    real t91, t92, t93, t94, t95, t96, t97, t98;
    real t99, t100, t101, t102, t103, t104, t105;
    real t106, t107, t108, t109, t110, t111, t112;
    real t113, t114, t115, t116, t117, t118, t119;
    real t120, t121, t122, t123, t124, t125, t126;
    real t127, t128;

    t1 = pow(6.0,0.333333333333333);
    t2 = 1/pow(3.141592653589793,0.333333333333333);
    t3 = 1/t1;
    t4 = 1/pow(3.141592653589793,2.666666666666667);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/pow(3.141592653589793,0.666666666666667);
    t8 = pow(rhoa,1.333333333333333);
    t9 = 1/t8;
    t10 = asinh(3.8978*t3*t7*grada*t9);
    t11 = 0.098225*t3*t7*grada*t10*t9;
    t12 = 4.166666666666667E-5*t3*t4*t5*t6+t11+1.0;
    t13 = 1/t12;
    t14 = pow(rhoa,0.333333333333333);
    t15 = 1/pow(6.0,0.666666666666667);
    t16 = 1/pow(3.141592653589793,1.333333333333333);
    t17 = pow(grada,2.0);
    t18 = 1/pow(rhoa,2.666666666666667);
    t19 = 1/pow(2.718281828459045,25.0*t15*t16*t17*t18);
    t20 = 0.2743-
        0.1508*t19;
    t21 = 0.25*t15*t16*t17*t18*t20+t11+1.0;
    t22 = 1/pow(rhoa,6.333333333333333);
    t23 = sqrt(15.19284484*t15*t16*t17*t18+1.0);
    t24 = 1/t23;
    t25 = 1/pow(rhoa,3.666666666666667);
    t26 = -0.510481873333333*t15*t16*t17*t24*t25;
    t27 = 1/pow(rhoa,2.333333333333333);
    t28 = -0.130966666666667*t3*t7*grada*t10*t27;
    t29 = t28+t26-2.222222222222222E-4*t3*t4*t5*t22;
    t30 = 1/pow(t12,2.0);
    t31 = t28+t26-0.666666666666667*t15*t16*t17*t20*t25-0.418888888888889*
        t3*t4*t5*t22*t19;
    t32 = pow(grada,3.0);
    t33 = 0.382861405*t15*t16*grada*t24*t18;
    t34 = 0.098225*t3*t7*t10*t9;
    t35 = t34+t33+1.6666666666666666E-4*t3*t4*t32*t6;
    t36 = 0.5*t15*t16*t18*t20*grada+t34+t33+0.314166666666667*
        t3*t4*t32*t6*t19;
    t37 = 1/pow(rhoa,0.666666666666667);
    t38 = pow(t29,2.0);
    t39 = 1/pow(t12,3.0);
    t40 = 1/pow(rhoa,7.333333333333333);
    t41 = 1/pow(t23,3.0);
    t42 = -1.723482643374637*t3*t4*t5*t41*t40;
    t43 = 1/pow(rhoa,4.666666666666667);
    t44 = 2.552409366666666*t15*t16*t17*t24*t43;
    t45 = 1/pow(rhoa,3.333333333333333);
    t46 = 0.305588888888889*t3*t7*grada*t10*t45;
    t47 = t46+t44+t42+0.001407407407407*t3*t4*t5*t40;
    t48 = 1/pow(3.141592653589793,4.0);
    t49 = pow(grada,6.0);
    t50 = 1/pow(rhoa,10.0);
    t51 = t46+t44+2.444444444444445*t15*t16*t17*t20*t43+t42-
        4.65432098765432*t48*t49*t50*t19+3.77*t3*t4*t5*t40*t19;
    t52 = 1.292611982530978*t3*t4*t32*t41*t22;
    t53 = -1.53144562*t15*t16*grada*t24*t25;
    t54 = -0.130966666666667*t3*t7*t10*t27;
    t55 = t54+t53+t52-8.888888888888888E-4*t3*t4*t32*t22;
    t56 = pow(grada,
        5.0);
    t57 = 1/pow(rhoa,9.0);
    t58 = -1.333333333333333*t15*t16*t20*t25*grada+t54+t53+
        t52+3.490740740740741*t48*t56*t57*t19-2.513333333333333*t3*
        t4*t32*t22*t19;
    t59 = pow(t35,2.0);
    t60 = -0.969458986898234*t3*t4*t17*t41*t6;
    t61 = 0.76572281*t15*t16*t24*t18;
    t62 = t61+t60+5.0E-4*t3*t4*t17*t6;
    t63 = 1/pow(rhoa,8.0);
    t64 = t61+t60+0.5*t15*t16*t18*t20-2.618055555555555*t48*
        t5*t63*t19+1.570833333333333*t3*t4*t17*t6*t19;
    t65 = 1/pow(rhoa,1.666666666666667);
    t66 = pow(t29,3.0);
    t67 = 1/pow(t12,4.0);
    t68 = 1/pow(t23,5.0);
    t69 = 1/pow(rhoa,11.0);
    t70 = -17.45640292348261*t48*t49*t68*t69;
    t71 = 1/pow(rhoa,8.333333333333334);
    t72 = 21.25628593495385*t3*t4*t5*t41*t71;
    t73 = 1/pow(rhoa,5.666666666666667);
    t74 = -13.49940953925926*t15*t16*t17*t24*t73;
    t75 = 1/pow(rhoa,4.333333333333333);
    t76 = -1.01862962962963*t3*t7*grada*t10*t75;
    t77 = t76+t74+t72-0.010320987654321*t3*t4*t5*t71+t70;
    t78 = 1/
        pow(3.141592653589793,5.333333333333333);
    t79 = pow(grada,8.0);
    t80 = 1/pow(rhoa,13.66666666666667);
    t81 = t76+t74-11.40740740740741*t15*t16*t17*t20*t73+t72+
        t70-310.2880658436214*t15*t78*t79*t80*t19-31.74246913580246*
        t3*t4*t5*t71*t19+88.43209876543209*t48*t49*t69*t19;
    t82 = 13.09230219261196*t48*t56*t68*t50;
    t83 = -13.35699048615344*t3*t4*t32*t41*t40;
    t84 = 6.295943104444444*t15*t16*grada*t24*t43;
    t85 = 0.305588888888889*t3*t7*t10*t45;
    t86 = t85+t84+t83+0.00562962962963*t3*t4*t32*t40+t82;
    t87 = pow(grada,
        7.0);
    t88 = 1/pow(rhoa,12.66666666666667);
    t89 = 4.88888888888889*t15*t16*t20*t43*grada+t85+t84+
        t83+t82+232.716049382716*t15*t78*t87*t88*t19-59.34259259259259*
        t48*t56*t50*t19+18.15185185185185*t3*t4*t32*t40*t19;
    t90 = -9.81922664445897*t48*t5*t68*t57;
    t91 = 7.755671895185867*t3*t4*t17*t41*t22;
    t92 = -2.041927493333334*t15*t16*t24*t25;
    t93 = t92+t91-0.002666666666667*t3*t4*t17*t22+t90;
    t94 = 1/pow(rhoa,11.66666666666667);
    t95 = t92+t91+t90-1.333333333333333*t15*t16*t20*t25-174.537037037037*
        t15*t78*t49*t94*t19+38.39814814814815*t48*t5*t57*t19-9.215555555555554*
        t3*t4*t17*t22*t19;
    t96 = pow(t35,3.0);
    t97 = 7.364419983344228*t48*t32*t68*t63;
    t98 = -3.877835947592934*t3*t4*grada*t41*t6;
    t99 = t98+0.001*t3*t4*grada*t6+t97;
    t100 = 1/pow(rhoa,10.66666666666667);
    t101 = 3.769999999999999*t3*t4*grada*t6*t19-23.5625*t48*
        t32*t63*t19+130.9027777777778*t15*t78*t56*t100*t19+t98+t97;
    t102 = 1/
        pow(t12,5.0);
    t103 = 1/pow(t23,7.0);
    t104 = 1/pow(rhoa,14.66666666666667);
    t105 = -1768.082807206625*t15*t78*t79*t103*t104;
    t106 = 1/pow(rhoa,12.0);
    t107 = 407.3160682145942*t48*t49*t68*t106;
    t108 = 1/pow(rhoa,9.333333333333334);
    t109 = -222.7122571383003*t3*t4*t5*t41*t108;
    t110 = 1/pow(rhoa,6.666666666666667);
    t111 = 81.79054014962962*t15*t16*t17*t24*t110;
    t112 = 4.414061728395062*t3*t7*grada*t10*t6;
    t113 = 1/pow(3.141592653589793,6.666666666666667);
    t114 = 1326.062105404968*t15*t78*t87*t103*t80;
    t115 = -266.2101445831098*t48*t56*t68*t69;
    t116 = 119.207549500079*t3*t4*t32*t41*t71;
    t117 = -30.96923364888889*t15*t16*grada*t24*t73;
    t118 = -1.01862962962963*t3*t7*t10*t75;
    t119 = -994.5465790537263*t15*t78*t49*t103*t88;
    t120 = 166.9268529558024*t48*t5*t68*t50;
    t121 = -56.0131859096757*t3*t4*t17*t41*t40;
    t122 = 7.487067475555555*t15*t16*t24*t43;
    t123 = 745.9099342902949*t15*t78*t56*t103*t94;
    t124 = -98.1922664445897*t48*t32*t68*t57;
    t125 = 20.68179172049565*t3*t4*grada*t41*t22;
    t126 = -559.4324507177213*t15*t78*t5*t103*t100;
    t127 = 51.5509398834096*t48*t17*t68*t63;
    t128 = -3.877835947592934*t3*t4*t41*t6;

   /* code */
    res[0] = 0.5*(-1.5*t1*t13*t2*t31*t8+1.5*t1*t2*t21*t29*
        t30*t8-2.0*t1*t13*t14*t2*t21);
    res[1] = 0.5*(1.5*t1*t2*t21*t30*t35*t8-1.5*t1*t13*t2*
        t36*t8);
    res[2] = 0.5*(-1.5*t1*t13*t2*t51*t8+1.5*t1*t2*t21*t30*
        t47*t8-3.0*t1*t2*t21*t38*t39*t8+3.0*t1*t2*t29*t30*t31*t8-0.666666666666667*
        t1*t13*t2*t21*t37-4.0*t1*t13*t14*t2*t31+4.0*t1*t14*t2*t21*
        t29*t30);
    res[3] = 0.5*(-1.5*t1*t13*t2*t58*t8+1.5*t1*t2*t21*t30*
        t55*t8-3.0*t1*t2*t21*t29*t35*t39*t8+1.5*t1*t2*t29*t30*t36*
        t8+1.5*t1*t2*t30*t31*t35*t8-2.0*t1*t13*t14*t2*t36+2.0*t1*t14*
        t2*t21*t30*t35);
    res[4] = 0.5*(-1.5*t1*t13*t2*t64*t8+1.5*t1*t2*t21*t30*
        t62*t8-3.0*t1*t2*t21*t39*t59*t8+3.0*t1*t2*t30*t35*t36*t8);
    res[5] = 0.5*(-1.5*t1*t13*t2*t8*t81+1.5*t1*t2*t21*t30*
        t77*t8+9.0*t1*t2*t21*t66*t67*t8+4.5*t1*t2*t29*t30*t51*t8-9.0*
        t1*t2*t21*t29*t39*t47*t8+4.5*t1*t2*t30*t31*t47*t8-9.0*t1*t2*
        t31*t38*t39*t8+0.444444444444444*t1*t13*t2*t21*t65-6.0*t1*
        t13*t14*t2*t51+6.0*t1*t14*t2*t21*t30*t47-12.0*t1*t14*t2*t21*
        t38*t39-2.0*t1*t13*t2*t31*t37+2.0*t1*t2*t21*t29*t30*t37+12.0*
        t1*t14*t2*t29*t30*t31);
    res[6] = 0.5*(-1.5*t1*t13*t2*t8*t89+1.5*t1*t2*t21*t30*
        t8*t86+9.0*t1*t2*t21*t35*t38*t67*t8+3.0*t1*t2*t29*t30*t58*
        t8-6.0*t1*t2*t21*t29*t39*t55*t8+3.0*t1*t2*t30*t31*t55*t8+1.5*
        t1*t2*t30*t35*t51*t8-3.0*t1*t2*t21*t35*t39*t47*t8+1.5*t1*t2*
        t30*t36*t47*t8-3.0*t1*t2*t36*t38*t39*t8-6.0*t1*t2*t29*t31*
        t35*t39*t8-4.0*t1*t13*t14*t2*t58+4.0*t1*t14*t2*t21*t30*t55-
        8.0*t1*t14*t2*t21*t29*t35*t39-0.666666666666667*t1*t13*t2*
        t36*t37+0.666666666666667*t1*t2*t21*t30*t35*t37+4.0*t1*t14*
        t2*t29*t30*t36+4.0*t1*t14*t2*t30*t31*t35);
    res[7] = 0.5*(-1.5*t1*t13*t2*t8*t95+1.5*t1*t2*t21*t30*
        t8*t93+9.0*t1*t2*t21*t29*t59*t67*t8+1.5*t1*t2*t29*t30*t64*
        t8-3.0*t1*t2*t21*t29*t39*t62*t8+1.5*t1*t2*t30*t31*t62*t8-3.0*
        t1*t2*t31*t39*t59*t8+3.0*t1*t2*t30*t35*t58*t8-6.0*t1*t2*t21*
        t35*t39*t55*t8+3.0*t1*t2*t30*t36*t55*t8-6.0*t1*t2*t29*t35*
        t36*t39*t8-2.0*t1*t13*t14*t2*t64+2.0*t1*t14*t2*t21*t30*t62-
        4.0*t1*t14*t2*t21*t39*t59+4.0*t1*t14*t2*t30*t35*t36);
    res[8] = 0.5*(1.5*t1*t2*t21*t30*t8*t99+9.0*t1*t2*t21*
        t67*t8*t96+4.5*t1*t2*t30*t35*t64*t8-9.0*t1*t2*t21*t35*t39*
        t62*t8+4.5*t1*t2*t30*t36*t62*t8-9.0*t1*t2*t36*t39*t59*t8-1.5*
        t1*t101*t13*t2*t8);
    res[9] = 0.5*(-1.5*t1*t13*t2*t8*(-3447.645176040237*t113*
        t19*t3*pow(grada,10.0)/pow(rhoa,17.33333333333333)+64.64197530864197*
        t110*t15*t16*t17*t20+283.6343209876543*t3*t4*t5*t108*t19-1325.447187928669*
        t48*t49*t106*t19+10136.0768175583*t15*t78*t79*t104*t19+t112+
        t111+t109+t107+t105)+6.0*t1*t2*t29*t30*t8*t81-8.0*t1*t13*t14*
        t2*t81-12.0*t1*t2*t21*t29*t39*t77*t8+6.0*t1*t2*t30*t31*t77*
        t8+36.0*t1*t2*t31*t66*t67*t8+54.0*t1*t2*t21*t38*t47*t67*t8+
        9.0*t1*t2*t30*t47*t51*t8-18.0*t1*t2*t38*t39*t51*t8-9.0*t1*
        t2*t21*t39*pow(t47,2.0)*t8-36.0*t1*t2*t29*t31*t39*t47*t8+1.5*
        t1*(t112+t111+t109+0.086008230452675*t3*t4*t5*t108+t107+t105)*
        t2*t21*t30*t8-36.0*t1*t102*t2*t21*pow(t29,4.0)*t8+8.0*t1*t14*
        t2*t21*t30*t77+48.0*t1*t14*t2*t21*t66*t67+1.777777777777778*
        t1*t13*t2*t31*t65-1.777777777777778*t1*t2*t21*t29*t30*t65-
        4.0*t1*t13*t2*t37*t51+24.0*t1*t14*t2*t29*t30*t51-48.0*t1*t14*
        t2*t21*t29*t39*t47+4.0*t1*t2*t21*t30*t37*t47+24.0*t1*t14*t2*
        t30*t31*t47-8.0*t1*t2*t21*t37*t38*t39-48.0*t1*t14*t2*t31*t38*
        t39+8.0*t1*t2*t29*t30*t31*t37-0.740740740740741*t1*t13*t18*
        t2*t21);
    res[10] = 0.5*(-1.5*t1*t13*t2*t8*(2585.733882030178*t113*
        t19*t3*pow(grada,9.0)/pow(rhoa,16.33333333333333)-22.81481481481481*
        t15*t16*t20*t73*grada-6903.909465020575*t15*t78*t87*t80*t19-
        141.3051851851852*t3*t4*t32*t71*t19+795.1131687242797*t48*
        t56*t69*t19+t118+t117+t116+t115+t114)+4.5*t1*t2*t29*t30*t8*
        t89-6.0*t1*t13*t14*t2*t89-9.0*t1*t2*t21*t29*t39*t8*t86+4.5*
        t1*t2*t30*t31*t8*t86+6.0*t1*t14*t2*t21*t30*t86+1.5*t1*t2*t30*
        t35*t8*t81-3.0*t1*t2*t21*t35*t39*t77*t8+1.5*t1*t2*t30*t36*
        t77*t8+9.0*t1*t2*t36*t66*t67*t8+27.0*t1*t2*t21*t38*t55*t67*
        t8+27.0*t1*t2*t21*t29*t35*t47*t67*t8+27.0*t1*t2*t31*t35*t38*
        t67*t8-36.0*t1*t102*t2*t21*t35*t66*t8+4.5*t1*t2*t30*t47*t58*
        t8-9.0*t1*t2*t38*t39*t58*t8+4.5*t1*t2*t30*t51*t55*t8-9.0*t1*
        t2*t21*t39*t47*t55*t8-18.0*t1*t2*t29*t31*t39*t55*t8-9.0*t1*
        t2*t29*t35*t39*t51*t8-9.0*t1*t2*t29*t36*t39*t47*t8-9.0*t1*
        t2*t31*t35*t39*t47*t8+1.5*t1*(t118+t117+t116-0.041283950617284*
        t3*t4*t32*t71+t115+t114)*t2*t21*t30*t8+36.0*t1*t14*t2*t21*
        t35*t38*t67+0.444444444444444*t1*t13*t2*t36*t65-0.444444444444444*
        t1*t2*t21*t30*t35*t65-2.0*t1*t13*t2*t37*t58+12.0*t1*t14*t2*
        t29*t30*t58-24.0*t1*t14*t2*t21*t29*t39*t55+2.0*t1*t2*t21*t30*
        t37*t55+12.0*t1*t14*t2*t30*t31*t55+6.0*t1*t14*t2*t30*t35*t51-
        12.0*t1*t14*t2*t21*t35*t39*t47+6.0*t1*t14*t2*t30*t36*t47-12.0*
        t1*t14*t2*t36*t38*t39-4.0*t1*t2*t21*t29*t35*t37*t39-24.0*t1*
        t14*t2*t29*t31*t35*t39+2.0*t1*t2*t29*t30*t36*t37+2.0*t1*t2*
        t30*t31*t35*t37);
    res[11] = 0.5*(-1.5*t1*t13*t2*t8*(-1939.300411522634*
        t113*t19*t3*t79/pow(rhoa,15.33333333333333)+4.88888888888889*
        t15*t16*t20*t43+4596.141975308642*t15*t78*t49*t88*t19-447.9783950617283*
        t48*t5*t50*t19+60.59925925925925*t3*t4*t17*t40*t19+t122+t121+
        t120+t119)+3.0*t1*t2*t29*t30*t8*t95-4.0*t1*t13*t14*t2*t95-
        6.0*t1*t2*t21*t29*t39*t8*t93+3.0*t1*t2*t30*t31*t8*t93+4.0*
        t1*t14*t2*t21*t30*t93+3.0*t1*t2*t30*t35*t8*t89-6.0*t1*t2*t21*
        t35*t39*t8*t86+3.0*t1*t2*t30*t36*t8*t86+9.0*t1*t2*t21*t38*
        t62*t67*t8+9.0*t1*t2*t21*t47*t59*t67*t8+18.0*t1*t2*t29*t31*
        t59*t67*t8+36.0*t1*t2*t21*t29*t35*t55*t67*t8+18.0*t1*t2*t35*
        t36*t38*t67*t8+1.5*t1*t2*t30*t47*t64*t8-3.0*t1*t2*t38*t39*
        t64*t8+1.5*t1*t2*t30*t51*t62*t8-3.0*t1*t2*t21*t39*t47*t62*
        t8-6.0*t1*t2*t29*t31*t39*t62*t8-3.0*t1*t2*t39*t51*t59*t8-36.0*
        t1*t102*t2*t21*t38*t59*t8+6.0*t1*t2*t30*t55*t58*t8-12.0*t1*
        t2*t29*t35*t39*t58*t8-6.0*t1*t2*t21*t39*pow(t55,2.0)*t8-12.0*
        t1*t2*t29*t36*t39*t55*t8-12.0*t1*t2*t31*t35*t39*t55*t8-6.0*
        t1*t2*t35*t36*t39*t47*t8+1.5*t1*(t122+t121+0.016888888888889*
        t3*t4*t17*t40+t120+t119)*t2*t21*t30*t8+24.0*t1*t14*t2*t21*
        t29*t59*t67-0.666666666666667*t1*t13*t2*t37*t64+4.0*t1*t14*
        t2*t29*t30*t64-8.0*t1*t14*t2*t21*t29*t39*t62+0.666666666666667*
        t1*t2*t21*t30*t37*t62+4.0*t1*t14*t2*t30*t31*t62-1.333333333333333*
        t1*t2*t21*t37*t39*t59-8.0*t1*t14*t2*t31*t39*t59+8.0*t1*t14*
        t2*t30*t35*t58-16.0*t1*t14*t2*t21*t35*t39*t55+8.0*t1*t14*t2*
        t30*t36*t55-16.0*t1*t14*t2*t29*t35*t36*t39+1.333333333333333*
        t1*t2*t30*t35*t36*t37);
    res[12] = 0.5*(-1.5*t1*t13*t2*t8*(1454.475308641976*t113*
        t19*t3*t87/pow(rhoa,14.33333333333333)-2967.12962962963*t15*
        t78*t56*t94*t19+230.3888888888889*t48*t32*t57*t19-20.10666666666666*
        t3*t4*grada*t22*t19+t125+t124+t123)-3.0*t1*t2*t21*t29*t39*
        t8*t99+1.5*t1*t2*t30*t31*t8*t99+2.0*t1*t14*t2*t21*t30*t99+
        9.0*t1*t2*t31*t67*t8*t96-36.0*t1*t102*t2*t21*t29*t8*t96+12.0*
        t1*t14*t2*t21*t67*t96+4.5*t1*t2*t30*t35*t8*t95-9.0*t1*t2*t21*
        t35*t39*t8*t93+4.5*t1*t2*t30*t36*t8*t93+27.0*t1*t2*t21*t29*
        t35*t62*t67*t8+27.0*t1*t2*t21*t55*t59*t67*t8+27.0*t1*t2*t29*
        t36*t59*t67*t8+4.5*t1*t2*t30*t55*t64*t8-9.0*t1*t2*t29*t35*
        t39*t64*t8+4.5*t1*t2*t30*t58*t62*t8-9.0*t1*t2*t21*t39*t55*
        t62*t8-9.0*t1*t2*t29*t36*t39*t62*t8-9.0*t1*t2*t31*t35*t39*
        t62*t8-9.0*t1*t2*t39*t58*t59*t8-18.0*t1*t2*t35*t36*t39*t55*
        t8+1.5*t1*t101*t2*t29*t30*t8+1.5*t1*(t125-0.005333333333333*
        t3*t4*grada*t22+t124+t123)*t2*t21*t30*t8+6.0*t1*t14*t2*t30*
        t35*t64-12.0*t1*t14*t2*t21*t35*t39*t62+6.0*t1*t14*t2*t30*t36*
        t62-12.0*t1*t14*t2*t36*t39*t59-2.0*t1*t101*t13*t14*t2);
    res[13] = 0.5*(-1.5*t1*t13*t2*t8*(-1090.856481481482*
        t113*t19*t3*t49/pow(rhoa,13.33333333333333)-102.1041666666666*
        t48*t17*t63*t19+3.769999999999999*t3*t4*t6*t19+1832.638888888889*
        t15*t78*t5*t100*t19+t128+t127+t126)-12.0*t1*t2*t21*t35*t39*
        t8*t99+6.0*t1*t2*t30*t36*t8*t99+36.0*t1*t2*t36*t67*t8*t96+
        54.0*t1*t2*t21*t59*t62*t67*t8+9.0*t1*t2*t30*t62*t64*t8-18.0*
        t1*t2*t39*t59*t64*t8-9.0*t1*t2*t21*t39*pow(t62,2.0)*t8-36.0*
        t1*t2*t35*t36*t39*t62*t8-36.0*t1*t102*t2*t21*pow(t35,4.0)*
        t8+6.0*t1*t101*t2*t30*t35*t8+1.5*t1*(t128+0.001*t3*t4*t6+t127+
        t126)*t2*t21*t30*t8);

}

static void
pw91x_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    pw91x_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
